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A  LOCALIZED  FINITE-ELEMENT  METHOD  FOR  THREE  DIMENSIONAL  SHIP  MOTION  PROBLEMS' 


Kwang  June  Bat 

David  W.  Taylor  Naval  Ship  Research  and  Development  Center 
Bethcada,  Maryland  200*4 


Abstract 

An  application  of  the  localized  finite-element  method 
to  a  three-dimensional  time-harmonic  free  surface  flow  in  a 
canal  is  presented.  Boundary  conditions  on  both  the  free 
surface  and  the  body  are  linearized  and  imposed  on  their 
equilibrium  positions.  By  utilizing  known  set  of  eigenfunc¬ 
tions,  the  computation  domain  is  reduced  to  a  very  small 
local  domain  where  an  eight -node  linear  three-dimensional 
dement  is  used.  Proper  matching  is  also  imposed  between 
two  sets  of  trial  functions  on  the  truncated  boundary.  To 
be  solved  are  the  problems  concerning:  (1)  six  degree-of- 
frtedom  radiation  and  diffraction  in  three  dimensions,  (2) 
two  dimensional  motion  corresponding  to  the  local  flow  at 
the  midship  cross-sectional  plane,  and  (J)  related  eigen¬ 
values.  Specifically,  two  sets  of  results  for  two  ship  loca¬ 
tions  in  a  canal  are  presented.  In  both  cases,  the  eigen¬ 
values  of  the  local  cross-sectional  plane  are  shown  to  play  a 
significant  role  in  the  three-dimensional  results.  A 
remarkable  similarity  between  exciting  forces  and  moment 
and  the  damping  coefficients  corresponding  to  the  modes 
of  their  motions  arc  also  observed  in  the  results.  The  ac¬ 
curacy  of  the  three-dimensional  results  presented  here  is 
also  discussed  by  comparing  two  sets  of  eigenvalues  com¬ 
puted  by  using  two  different  sizes  of  finite  elements. 

1.  Introduction 

Steady-state  time-harmonic  motions  of  an  invisetd,  In¬ 
compressible  fluid  with  free  surface  in  the  presence  of  a 
body  or  bodies  in  it  arc  described  by  a  boundary-value  pro¬ 
blem  governed  by  the  Laplace  equation  with  appropriate 
boundary  conditions.  In  the  past,  problems  of  this  type 
were  generally  solved  by  distributing  sources  (and/or 
dipoles)  on  the  body  boundary  and  using  Green's  theorem 
to  obtain  an  integral  equation  for  the  strength  of  these 
boundary  singularities:  or,  alternatively,  by  using  sources 
and  higher-order  multipole  expansions  at  an  interior  point 
within  the  body,  the  strengths  of  these  singularities  being 
determined  so  as  to  satisfy  the  body  boundary  condition. 

In  all  cases  it  is  conventional  to  utilize  the  singularities 
which  are  solutions  of  the  boundary-value  problem  stated 
above,  except  for  the  body  boundary  condition  which  is  in¬ 
voked  separately  to  determine  the  singularity  distribution. 


This  work  was  supported  by  the  Numerical  Naval  Hydro¬ 
dynamics  Program  at  the  David  W.  Taylor  Naval  Ship 
R&D  Center.  This  Program  is  jointly  supported  by 
NSRDC  and  the  Office  of  Naval  Research.  Part  of  the 
present  work  was  done  while  the  author  was  a  visiting 
professor  at  the  Ecole  Nationale  Superieure  dc  Techniques 
Avancecs  (ENSTA)  in  1980.  and  partially  supported  by 
the  Ministry  of  Defense,  France. 


For  two-  and  three-dimensional  motions  in  a  fluid  of  in¬ 
finite  depth,  or  of  finite  but  constant  depth,  the  required 
singularities  are  well  known,  although  of  rather  com¬ 
plicated  analytical  form,  so  that  the  approach  described  in 
the  foregoing  corresponds  to  solving  a  Fredholm  integral 
equation  over  the  body  surface,  with  a  rather  complicated 
kernel  function.  An  extensive  list  of  literature  on  this  sub¬ 
ject  can  be  found  in  Wehausen1. 

In  this  paper  a  numerically  oriented  as  well  as  more 
versatile  method  is  introduced  as  an  alternative  approach  to 
the  solution  of  the  problem.  This  alternative  approach  is 
based  on  a  variational  principle  which  is  utilized  to  deter¬ 
mine  the  velocity  potential  throughout  the  fluid  domain. 
The  present  work  is  a  direct  extension  of  the  earlier  work 
by  Bai  and  Yeung2  for  a  general  three-dimensional  body 
geometry  in  a  canal.  In  this  procedure,  the  known  solution 
space  in  certain  subdomains  is  made  use  of  in  order  to 
reduce  the  ‘computation  box*  to  a  size  as  small  as  possible. 
We  call  this  drastically-reduced  computation  box  the 
localized  finite-element  domain.  The  phrase  “localized 
finite-clement  method"  is  used  to  denote  the  finite  element 
method  applied  only  in  this  localized  subdomain.  In  this 
particular  problem,  i.e.,  in  a  free-surface  flow  in  a  canal, 
the  solution  space  is  represented  by  the  complete  set  of 
eigenfunctions  in  the  subdomain.  It  is  essential  to  make  the 
domain  of  computation  as  small  as  possible  because  the 
practical  success  or  failure  of  a  numerical  method,  in 
general,  mainly  depends  upon  the  size  of  the  computation 
box.  In  the  present  procedure,  the  reduction  of  the  original 
infinite  fluid  domain  to  a  small  localized  subdomain  is 
achieved  by  replacing  the  conventional  radiation  condition 
by  the  matching  condition.  This  requires  that  the  velocity 
potential  and  its  normal  derivative  (i.e„  normal  velocity) 
represented  by  two  sets  of  different  trial  functions,  defined 
in  the  adjacent  subdomains,  be  continuous  along  the  fic¬ 
titious  interface  juncture  boundary  surface.  Then  we  ap¬ 
proximate  the  potential  in  (he  localized  finite-element  do¬ 
main  by  piecewise  polynomial  trial  functions  defined  in 
each  finite  element.  Thus,  in  effect,  an  integral  equation 
over  a  body  surface  with  a  complicated  kernel  is  replaced 
by  a  system  of  equations  over  a  much  larger  fluid  domain, 
but  with  a  much  simpler  kernel. 

Previous  studies  of  the  finite  element  method  applied 
to  time-harmonic  ship  motion  problems  have  been  made  by 
Bai34,  Berkhoff7,  Smith*,  Chen  and  Mei9,  Seto  and 
Yamamoto10,  Yue  et  al11,  Euvrard  cl  aJ*2,  and 
Chowdbury13. 

Additional  numerical  results  of  the  related  two- 
dimensional  boundary-value  and  eigenvalue  problems  at  the 
midship  cross-sectional  plane  are  also  presented.  Two  dif¬ 
ferent  sizes  of  the  finite  elements  in  the  midship  cross- 
sectional  plane  are  used  in  these  additional  computations. 
These  results  arc  used  as  reference  in  discussing  (he  ac¬ 
curacy  of  the  present  three-dimensional  results. 
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H.  Mathematical  Formulation 


with 


Considered  here  is  steady-state  time-harmonic  f ree- 
surface  flow  in  the  presence  of  a  body  floating  or  submerg¬ 
ed  in  a  canal  with  a  rectangular  uniform  cross-section. 
However,  a  local  variation  of  the  bottom  and  side-wall 
geometry  is  present  and  a  uniform  rectangular  cross-section 
of  each  side  (not  necessarily  the  same  on  both  sides),  can 
be  similarly  treated  by  the  present  method.  The  coordinate 
system  is  right-handed  and  rectangular ,  The  y-axis  is 
directed  opposite  to  the  force  of  gravity,  and  the  xz-plane 
coincides  with  the  undisturbed  free  surface.  The  bottom  of 
the  canal  is  in  the  y  «  -H  plane  and  the  side  walls  are  in 
the  2  ~  x  W/2  planes.  We  neglect  surface  tension,  and 
assume  that  the  fluid  is  inviscid  and  incompressible  and 
that  the  motion  is  irrotational.  Furthermore,  we  assure  that 
the  motion  of  the  body  and.  consequently,  the  generated 
wave,  to  be  small  in  some  sense,  so  that  the  boundary  con¬ 
dition  on  the  body  and  on  the  free  surface  can  be  lineariz¬ 
ed  and  satisfied  at  the  mean  equilibrium  positions  instead 
of  at  their  instantaneous  positions. 


(nlt  n2,  nj)  -  (n„.  n^,  n*) 

(n*  nJt  Dg)  »  (T-7j  a. 

where  IT  is  the  unit  normal  vector  into  the  body  with  com¬ 
ponents  (nx,  riy,  nz),  rc  is  the  vector  from  the  origin  to  the 
center  of  rotation,  and 


For  m  =  7, 


r  -  (*,  y.  a). 


(7) 


where  it  is  the  spatial  potential  associated  with  the  incident 
wave  system.  The  incident  wave  potential  of  unit  amplitude 
incoming  frorA  x  -  -<*>  is  given  by 


♦|  (*.  y.  x)  - 


g  cosh  m„  (y  +  H) 
0  cosh  tn0  H 


<») 


Let  ♦  (x,  y,  z.  t)  be  the  velocity  potential  describing 
the  flow  field.  The  continuity  equation  requires  <t>  to  satisfy 
Laplace’s  equation.  Let  o  be  the  angular  frequency  of  the 
time-harmonic  solution.  Then,  introducing  the  usual  lime 
and  spatial  decomposition  we  have 


where  m0  is  wave  number.  Finally,  to  make  the  solution  of 
this  problem  unique,  we  impose  the  radiation  condition  re¬ 
quiring  that  the  generated  waves  must  be  outgoing. 

III.  Localized  Finite  Element  Method 


♦  (x,  y ,  z.  t)  -  rteO<m>  (x,  y,  z)  S<m>(t)]  (1) 


mm  I,  ..  ,7 


where  ♦<"*>  -  ♦<m)  +  i*<m)  (2) 

is  the  complex-valued  spatial  potential,  and  o<'">(t)  is  the 
complex  time-harmonic  motion  amplitude  of  the  body  cor¬ 
responding  to  the  m-th  mode,  i.e.. 


•m(t)  »  lfe[«<'")(t)] 

-  Re  +  to<,m))e-'<«j  & 

where  em(t)  is  the  body  motion  amplitude.  For  m  *  7,  the 
diffraction  problem  corresponding  to  an  incoming  wave 
system  of  unit  amplitude,  one  simply  sets  a)7'  =  0  and  oj,7» 
■»  l/o  in  Equation  (3).  The  potential  function  +  (x,  y,  z) 
must  satisfy 


(&  +  m  +  B)*(x'y'z) 


in  the  fluid. 


<«> 


The  fluid  domain  defined  in  the  foregoing  formulation 
is  unbounded  along  the  x-axis.  For  numerical  computa¬ 
tions,  it  is  highly  desirable  that  the  computation  box  be 
made  as  small  as  possible.  The  goal  of  reducing  the 
original  infinite  fluid  domain  to  a  manageable  finite  do¬ 
main  is  achieved  by  making  use  of  the  known  solution 
space  in -the  truncated  infinite  subdomains  which  will  be 
defined  later.  As  a  result,  the  computation  domain  is 
reduced  to  a  very  local  subdomain,  called  a  localized  finite 
element  domain,  which  may  barely  include  any  source  of 
disturbance  in  the  fluid.  The  present  numerical  method  is 
called  a  localized  finite  clement  method  because  the  finite 
element  numerical  computations  are  made  only  for  a  local 
domain.  Throughout  this  section,  the  subscript  m  in  the 
potential,  defined  in  the  Mathematical  Formulation  section 
is  omitted  with  the  understanding  that  +  will  always  depend 
upon  m. 

Let  us  draw  two  imaginary  vertical  planes  J(  and  J2 
which  separate  the  original  fluid  domain  D  into  the  three 
subdomains:  D0,  D(,  D2  as  shown  in  Figure  1.  We  assure 
that  D0  includes  the  ship  (and/or  any  other  sources  of 
disturbance).  The  boundary  surfaces  dD0,  dDj,  and  3D2 
are  denoted,  respectively,  as 


V*lsF-c 

Mso  -  v«  -  f(5> 


(5a) 

<3b) 


9D0  m  Sfg  +  Swo  +  Sgo  +  j|  +  jg  +  So 
dDj  m  5p  +  Syyj  +  Sg  +  Jj  +  Sjy  i  *  1,2 


(9) 


<*> 

Ms."0  <*> 

where  v  «  o*/g 
g  —  acceleration  due  to  gravity 

SF  —  undisturbed  free  surface,  y  ■  0,  outside  of  the  body 
Sb  —  bottom  surface 
S*  -  side  walls 

So  —  body  surface  below  the  mean  free  surface. 

The  normal  velocity  f(s)  depends  upon  the  mode  index  m. 

For  m  »  1,2 . 6.  corresponding  to  the  sway,  heave, 

surge,  roll,  pitch,  and  yaw  modes  of  motion,  respectively, 
f(s)  is  given  by 


where  Sp,  Sgj,  and  Swj.  denote,  respectively,  the  free  sur¬ 
face,  the  bottom,  and  the  canal  side  walls,  in  the  subdo¬ 
mains  Dj  (for  i  •0,1,2)  and  where  Ssj  (for  i  »  1,  2)  are 
the  boundaries  at  infinity.  The  ship-hullsurface  is  denoted 
by  V 


Let  +0,  4|.  and  +2  denote  the  velocity  potentials  defined 
in  the  subdomains  D0,  D(,  and  D2,  respectively.  Then  we 
have,  from  Equations  (4)  and  (5), 

V^p  -  0  in  D„ 


♦on  -  *♦*  -  0  on  Sfp 
♦on  -  «»)  on  Sp 


(10) 


ff')(s)  -  nj 


(6) 


♦on  “0 


°n  $BsU  Swo 


2 


1 


Figure  1  —  Boundary  Configurations  of  the  Three  Subdivided  Fluid  Domains 


and,  for  4],  i  -  1,  2, 

V2  f,  -  0  in  D, 

4i„-i4j-0  on  SR  (II) 

♦in  ”  0  on  Sgjij  S^fj 

It  is  understood  in  Equation  (11)  that  a  proper  radiation 
condition  is  imposed  at  infinities  for  |j  and  In  addition, 
we  require 


KlWoil.+l)  *  //l  jfr+o)2*'-  %,  ♦o2* 

-ffs  f(i#o*  »o- 1  Mm*  04> 

*JJl  (*o-  jMin* 


♦o  ”  4i 

♦on  +  ♦in  "  0 


on  Jj,  i  —  1.  2 


(12) 


where  the  normal  vector  is  taken  outwards  from  the  fluid 
subdomain  defined  for  each  potential,  e.g.. 


♦la  ”  ♦!*•  +00  “  -4o*  on  J1 

+2n  *  -42*.  +00  »  +0*  on  h 


(13) 


and 


K2(M|.+2l  “  fffD  to7** 

M,‘4oHon+  2 


By  the  juncture  conditions  in  Equations  (12),  the  solutions 
of  Equations  (10)  through  (12)  are  unique  and  identical  to 
the  solution  of  the  original  formulation  in  the  Mathematical 
Formulation  section.  (This  can  be  shown  by  applying 
Green's  theorem  to  the  difference  of  the  solution  of  the 
original  equation,  defined  in  the  previous  section,  and  the 
solutions  ♦„,  ♦),  and  +2.  defined  in  the  three  subdomains). 

Let  us  assume  that  the  general  solutions  of  Equation 
(II)  with  appropriate  radiation  conditions,  are  known,  i.e„ 
the  compiete  set  of  eigenfunctions  or  Green's  functions  of 
the  problems.  Then  the  first  step  toward  the  construction 
of  the  appropriate  functional  for  our  variational  equation 
equivalent  to  the  coupled  partial  differential  equations 
defined  in  (10)  through  (12)  is  the  choice  of  trial  function 
space  for  each  function  >0,  and  +2-  If  the  bases  of  the 
trial  functions  for  4t  and  |2  are  chosen  from  the  solution 
space  (the  known  general  solutions),  then  we  can  obtain 
the  following  two  functionals: 


+jfjT  ((♦r^oMon+yMlnl*  <l5) 


Setting  the  tint  variation  of  either  functional  to  be  zero,  i.e., 
♦K,  {♦„.♦,.  42}  -  0  (16a) 

or  *  *2  <♦»•  4|.  -  0  (16b) 

is  the  same  as  solving  Equations  (10)  through  (12).  It  is 
easy  to  show  this  by  taking  the  first  variation  and  by  using 
Green's  theorem2.  The  integral  expressions  of  the  func¬ 
tionals  involve  only  the  subdomain  D0  which  we  call  the 
“localized  finite-element  domain.”  If  one  takes  the  localiz¬ 
ed  finite-element  domain  to  be  small,  then  the  domain  over 
which  the  integrals  have  to  be  computed  will  also  be  small. 
On  the  other  hand,  one  has  to  take  many  terms  (eigenfunc- 


3 


lions)  to  represent  the  trial  solutions  4|  and  +2  in  the  com¬ 
putation  of  the  approximate  solutions,  and  vice  versa. 


reduces  to  a  set  of  linear  algebraic  equations.  (In  this  pro¬ 
cedure,  only  the  coefficients  are  subject  to  variation.) 


It  is  interesting  that  the  stationary  value  of  both  func¬ 
tions  are 

K|{+o.*|.*2> 

T  f(s)+ods 
'S0 

(17a) 

*2  *2)  --W 

(17b) 

i>  *0,  ♦(.  and  are  the  exact  solutions.  This  stationary 
value  is  simply  <-l/2e)  times  the  added  mass  and  (-l/2qo) 
times  the  damping  coefficients. 


Either  functional,  defined  in  (14)  or  (15),  will  give  a 
variational  equation  mathematically  equivalent  to  the 
Original  boundary -value  problem  at  hand.  The  use  of  the 
functional  given  in  (14)  is  slightly  advantageous  because  the 
norma)  derivative  of  ♦„  is  not  involved  in  the  coupling  in¬ 
tegral  terms.  In  the  present  computation,  the  functional 
given  in  Equation  (14)  is  used.  The  equivalence  of  the  dif¬ 
ferential  equation  to  the  variational  problem  is  basic  to  the 
choice  of  the  computational  scheme.  One  significant  dif¬ 
ference  between  the  functional  method  and  the  differentia] 
equation  is  the  fact  that  the  expressions  for  the  associated 
functionals  in  (14)  and  (15)  involve  no  second  derivatives, 
owing  to  the  integration  by  parts  used  to  construct  these 
functionals,  i.c..  Green’s  theorem  is  used  here.  It  follows 
(hat  the  functionals  will  be  well  defined  if  only  the  first 
derivative  of  the  function,  rather  than  the  second,  is  re¬ 
quired  to  be  bounded.  Therefore,  the  class  of  admissible 
functions,  in  the  problem  to  find  the  stationary  point,  is 
enlarged  to  a  space  bigger  than  that  for  the  original  dif¬ 
ferential  equation.  We  now  have  the  advantage,  while  sear¬ 
ching  for  the  stationary  point  of  the  functional,  of  being 
permitted  to  try  functions  outside  the  class  of  those 
originally  admissible.  In  practice,  this  means  that  we  can 
now  try  continuous  functions  whose  first  derivatives  are 
only  piecewise  continuous.  In  other  words,  the  first 
derivative  can  have  finite  discontinuities  at  the  juncture 
boundary  between  adjacent  elements.  It  is  very  easy  to  con¬ 
struct  basis  functions  that  satisfy  the  previous  re¬ 
quirements. 

In  the  localised  finite  element  domain  D0,  the  basis  for 
the  trial  function  is  chosen  from  a  polynomial  basis. 
Specifically,  eight-node  isoparametric,  linear  three- 
dimensional  elements  were  used  in  the  present  numerical 
calculations.  As  mentioned  earlier,  the  trial  functions  in  D| 
and  Dj  will  be  chosen  from  a  subspace  of  the  solution 
space  which  satisfies  the  Laplace  equation  with  the  free- 
surface  condition,  the  side-wall  condition,  the  bottom  con¬ 
dition,  and  the  radiation  condition  at  infinity.  The  eigen¬ 
functions  or  the  Green  functions  of  the  above  problem  can 
represent  the  solution  space.  However,  we  will  choose  the 
eigen-function  space  in  this  paper  for  its  simplicity. 

The  variational  equation  (16a)  goes  into  operational 
form  in  the  following  way.  Let  qg  (i  -  0,  1,  2,  and  j  - 
1,2,...  M^)  be  the  basis  for  the  trial  functions  in  each 
subdomain  Oj  (i  »  0,  I,  2).  Then  the  solution  is  assumed 
lobe 
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in  Dj  ft  »  0,  I.  2),  where  «pg  art  coefficients  to  be  deter¬ 
mine!.  By  substituting  Equation  (It)  in  the  functional 
defined  in  Equation  (14),  the  variational  equation  (16a) 


A  complete  set  of  eigenfunctions  (the  resonance  fre¬ 
quencies  are  excluded)  are  given  in  Wehausen14  as 


ITi  KdqX  cosh  m0(y  4  H)  nn  /  W\ 
r  cosh  m0H  W  V  2/' 

e*  cos  mp  (y  +  H)  cos  (z  -  yJJ 


(19) 


where  W  is  the  width  of  the  tapk  and  the  upper  and  lower 
signs  are  to  be  taken  in  the  subdomains  Dj  and  Dj,  respec¬ 
tively,  and  where ' 
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‘np  -  [V  -  (m 

and  m0  and  mp  ( p  —  I,  2,  .  .  .)  are  the  real  roots  of 
mg  tanh  m^H  **  v 
nip  tan  mpH  -  -v 

The  exponent  of  the  first  term  in  Equation  (19)  becomes 
real  when  nn/W  >01,,,  resulting  in  a  local  disturbance,  and 
if  n  -  0,  it  becomes  purely  a  two-dimensional  case.  The 
case  when  mQ  =  nn/W  (n  *  1,2,...),  i.e.,  the  case  of 
resonance,  is  left  out  in  the  present  study.  Here  we  con¬ 
sidered  only  the  case  of  m„  *  nn/W. 

IV.  Results  and  Discussions 

After  the  potential  +(■>,  defined  in  the  Mathematical 
Formulation  section,  has  been  computed,  the 
hydrodynamic  coefficients  can  be  computed  by 


♦°>nids 


(22) 


where  jig  are  the  added  masses  for  i,  j  =  1,  2,  3,  moments 
of  added  mass  for  i  (or  j)  =  1,  2,  3  and  j  (or  i)  =  4,  5,  6, 
and  moments  of  inertia  of  added  mass  for  i,  j  •  4,  5.  6, 
and  where  Ag  (i,  j  -  I.  ...  6)  are  the  damping  coeffi¬ 
cients.  in  presenting  our  results,  the  nondimensional 
hydrodynamic  coefficients,  jig  and  Ag  are  defined  by  using 
the  ship  displacement  V  and  the  ship  length  L  as  follows: 


Pjj  -  <V«V’ 

ig  -  VeoV 


j  -  I,  2.  3 


Ftj  “  Pg/eVL,  i  -  1.2.3.  j -4.3.6 

«r  (23) 

Ag  -  Ig/goVL  i  -  4.  5.  6,  j  -  I,  2,  3 


Mij  -  Pg/eVL1, 

i.j-4.3,6 

kg  -  Ig/poVL* 


where  V  -  L  B  T  is  the  ship  displacement  volume.  The 
results  of  wave  excitation  forces  and  moments  are  also 
similarly  nondimensionalized  by  using  the  unit  amplitude 


4. 


i 


of  the  incoming  wave,  Y  =  1,  and  the  ship  length  L  as 


Fj  -  Fj/ffgYL*  ) 

_  Vi  =  I.  2.  3 

Mj  -  Mj/pgYL5  ) 


(24) 


where 
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In  the  present  computations,  a  rectangular  barge  which 
floats  parallel  to  the  tank  walls,  i.e.,  the  angle  between  the 
center  plane  of  the  barge  and  the  xy -plane  is  zero,  is 
treated  for  its  simplicity  in  data  preparation.  However,  any 
arbitrary  angle  and  ship  geometry  can  be  handled  by  the 
present  computer  program.  The  midship  section  plane  is 
assumed  to  coincide  with  the  yz-plane  for  the  present  test 
model.  The  cross-section  of  the  canal  at  the  origin,  i.e.,  the 
yz-plane  is  given  in  Figure  2.  We  denote  the  distance  bet¬ 
ween  the  ship  center-plane  and  the  x>  -plane  by  s,  as  shown 
in  Figure  2.  T  he  moment  is  computed  with  respect  to  the 
centroid  of  the  ship’s  water-plane  area.  i.e. ,7^  =  (0,  0,  s) 
in  Equation  (6).  In  presenting  our  numerical  results,  the 
nondimension al  wave  number,  vB  =  o-B/g  is  plotted  along 
the  abscissa. 


In  the  present  computations,  two  locations  of  the  body 
in  a  canal  are  computed;  first,  the  ship  is  located  at  the 
center  of  the  canal,  i.e.,  s/B  =  0,  and  second,  the  ship  is 
at  an  off-center  location  in  a  canal.  In  both  cases,  the  ship 
is  parallel  to  the  canal  walls,  as  mentioned  earlier. 
Specifically,  the  computations  are  made  for  W/B  *  H/T 
•  W/H  •=  B/T  =  2  and  LT  =  10  for  both  values  of 
s/B  «  0  and  0.125.  Before  we  compute  these  two  cases,  we 
have  done  computations  for  a  two-dimensional  rectangular 
cylinder  which  is  uniform  between  two  walls,  i.e.,  the 
cylinder  lies  normal  to  the  xy-plane.  These  test  results  were 
compared  with  those  computed  by  the  results  of  our 
previous  work  in  two-dimensional  problems.  The  com- 
oarison  was  good. 

Several  sets  of  finite  element  subdivisions  were  initially 
tested  in  the  present  study.  The  total  number  of  nodes 
tested  are  180,  350,  507,  1404,  and  1782  (the  corresponding 
numbers  of  nodes  on  the  snip  boundary  are,  44,  81,  81, 

249,  and  249,  respectively).  The  indices  n  in  the  eigenfunc¬ 
tion  (Equation  (19))  were  taken  between  6  and  10,  and 
while  the  indices  p  were  taken  between  6  and  10  for  n  =  0, 
the  values  of  p  were  reduced  as  n  increased  The  present 
computations  were  made  for  the  range  of  non -dimensional 
wave  numbers,  2n/m0L  -  0.52  through  2.72.  This  range 


was  selected  partly  due  to  the  limitation  of  numerical  ac¬ 
curacy  without  taking  very  small  finite  elements,  while  the 
first  mode  of  canal  resonance,  i.e.,  when  m0  W/»  -  I  (see 
the  text  following  Equation  (21)),  is  included.  This  first 
mode  of  the  canal  resonance,  mD  W/x  •  I,  corresponds  to 
vB  -  i. 44066  in  the  present  ship-canal  geometry. 

Considering  the  wave  number  range  to  be  tested  and 
the  computation  time,  we  decided  to  use  the  307-node 
model  throughout  the  present  computations.  In  subdividing 
the  localized  rmite-elcment  domain  into  a  number  of  finite 
elements,  we  used  four  equal  elements  along  the  depth  and 
four  equal  elements  across  the  width  of  the  canal.  Five 
equal  elements  were  used  along  both  sides  of  the  x-axis  and 
five  more  elements  were  used  further  along  both  sides  of 
the  x-axis.  The  ship  replaces  the  middle  ten  elements  along 
the  x-axis  and  the  middle  two  dements  along  the  z-axis  and 
the  two  elements  along  the  depth  from  the  free  surface. 
When  the  ship  is  at  an  off-center  location  in  a  canal,  we 
simply  stretch  the  elements  on  one  side  and  shrink  those  on 
the  other  side  mainly  for  the  ease  of  data  preparation, 
especially  for  the  case  when  the  off-center  distance  is  not 
too  great.  Any  complicated  situation  can  be  accommodated 
by  the  present  method  but  only  through  more  tedius  data 
preparation. 

For  all  cases  treated  here,  the  computations  were  made 
for  23  values  of  vB  between  0.2  and  2.4  with  a  constant  in¬ 
crement  of  0. 1  throughout  the  present  computations.  Then 
the  23  computed  values  of  the  hydrodynamic  coefficients 
were  plotted  by  a  computer  using  a  straight-line 
interpolation. 

The  computational  results  for  the  case  when  the  ship  is 
located  at  the  center  of  a  canal  are  presented  in  Figures  3 
through  6.  Figure  7  shows  the  hydrodynamic  coefficients 
computed  for  the  midship  section,  i.e.,  the  yz-plane  as  a 
purely  two  dimensional  problem  with  no  wave  radiation 
present.  The  result  for  a  two-dimensional  model  provides 
the  resonance  mode  which  plays  a  significant  role  in  the 
local  (low  field  in  three-dimensional  problem. 

In  Figure  3,  the  computed  values  of  added  mass  and 
moment  of  inertia  of  added  mass  are  shown.  The  surge 
added-mass  coefficient,  p j g,  shows  an  oscillatory  behavior, 
being  approximately  zero  at  vB  «  0.9  and  reaching  a  local 
maximum  at  vB  «  1.4.  The  heave  added-mass  coefficient 
M22  is  negative  approximately  between  vB  «=  0.20  and  1.20. 
fiioth  sway  added-mass  coefficients  pjj  and  roll  moment  of 
inertia  of  added  mass  P44  change  abruptly  from  positive 
peak  values  to  negative  peak  values  between  vB  -  0.9  and 
1.0.  A  similar  behavior  is  also  observed  in  the  yaw  moment 
of  inertia  of  added  mass  pjj  between  vB  *■  1.2  and  1.3. 

The  pitch  moment  of  inertia  of  added  mass  pgg  is 
minimum  around  vB  =  1.3. 

In  Figure  4,  the  six  diagonal  damping  coefficients,  ky 
(i  *=  1,  6),  are  shown  for  s/B  -  0.  The  value  of  1||  de¬ 
creases  with  increasing  value  of  vB.  and  is  approximately 
zero  at  vB  «  1.2  and  reaches  a  local  maximum  at  1 .8.  The 
values  of  >22  decrease  monotonically  and  also  rapidly  to 
zero  at  vB  »  I.  The  values  of  kjj  and  A44  are  zero  between 
vB  «  0.2  and  1.4,  then  sharply  increase  to  local  maxima  at 
2.0,  and  finally  decrease.  The  value  of  kjj  shows  a  behavior 
similar  to  kjj  and  144.  but  becomes  zero  again  at  vB  «  2.0 
and  increases  again.  The  pitch  damping  coefficient  kgg de¬ 
creases  rather  rapidly  to  zero  between  vB  »  1.0  and  1.8.  It 
is  of  interest  to  note  that  the  values  of  kjj.  I44,  and  kjj  are 
zero  for  vBS  1.4.  This  is  because  the  (low  fields  for  the 
sway,  roll,  and  yaw  motions  are  purely  due  to  local  distur¬ 
bances  when  vB<  */2  tanh  n/2  -  1.4406  -  -(this  value 
corresponds  to  the  first  mode  over  the  canal  resonance). 
(Note  that  the  motions  are  asymmetric  with  respect  to  the 
xy-plane  while  the  only  admissible  far -Field  wave  solution  is 
pure  two-dimensional). 
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Figure  3  —  Added  Mums,  and  Moments  of  Inertia  of  Added  Mass,  ^  emus  vB,  for  s/B  ■=  0 


Due  to  the  symmetry  of  the  body  and  the  canal 
geometry  with  respect  to  the  xy-plane  and  the  yz-plane, 
there  are  only  four  nonzero  coupling  hydrodynamic  coeffi¬ 
cients.  These  are  shown  in  Figure  5.  The  pitch  moment  of 
surge  added  mass  pfg  reaches  a  maximum  at  vB  «  1.1  and 
then  decreases.  The  roll-sway  coupling  term  h]4  changes 
abruptly  from  a  negative  peak  to  a  positive  peak  between 
»B  «  0.9  and  1.0.  The  surge-pitch  coupling  damping  coef¬ 
ficient  1  |g  is  a  maximum  at  vB  -  1.4  and  slowly  decreases 
to  a  positive  constant.  The  sway -roll  coupling  damping  coef¬ 
ficient  lj4  is  zero,  u  discussed  earlier,  between  vB  =  0.2 
and  1.4,  reaches  a  negative  peak  sharply,  and  then  increases. 

In  Figure  4,  the  magnitudes  of  the  wave  excitation 
forces  and  moment  are  shown.  The  magnitude  of  the  surge 
excitation  force  |F||  reaches  a  local  minimum  at  vB  =  1.2 
and  a  local  maximum  at  1.6.  The  magnitude  of  the  heave 


excitation  force  |F2|  decreases  monotonically  from  its  value 
at  vB  =  0.2  to  zero  at  1.2  and  then  reaches  a  local  max¬ 
imum  at  1.8.  It  is  of  interest  to  note  that  the  zero  value  of 
the  heave  excitation  force  would  be  at  vB  *=  1.0683,  which 
gives  an  incident  wavelength  equal  to  the  ship  length,  if  the 
Froude-Krylov  approximation  is  made;  this  value  is  not  too 
far  from  the  computed  result  of  vB  =  1.2.  The  magnitude 
of  the  pilch  excitation  moment  decreases  from  a  local  max¬ 
imum  of  zero  at  vB  =  1.8  and  slowly  increases  to  the  next 
local  maximum.  The  rest  of  the  other  excitation  forces  and 
moments,  not  shown  here,  are  all  zero,  due  to  symmetry  in 
the  problem. 

It  should  be  noted  in  comparing  the  damping  coeffi¬ 
cients  li|,  122-  and  *66  in  Figure  4,  the  excitation  forces 
|F||  and  |F2|,  and  the  moment  |Mj|,  in  Figure  6,  respec¬ 
tively.  There  exists  a  close  similarity  in  the  behavior  of  the 
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Figure  4  —  Damping  Coefficients,  ig  emus  vB,  for  s/B  -  ft 


damping  coefficients  and  that  of  wave  excitations  cor¬ 
responding  to  each  mode  of  motion  in  the  range  of  the 
values  vB  =  0.2  through  1.44  This  is  not  a  coincidence 
because  there  exist,  for  vB  <  I  44,  only  a  purely  two- 
dimensional  radiating  wave  (see  the  text  following  Equation 
(21)).  On  the  other  hand,  if  the  wavelength  becomes  short 
and  if  there  exist  many  modes  of  the  three-dimensional 
wave  system,  then  this  similarity  may  no  longer  hold. 

In  Figure  7,  the  hydrodynamic  coefficients  computed 
from  a  purely  two-dimensional  problem  in  the  yz-plane  (x 
■  0)  are  presented  as  mentioned  earlier.  These  computa¬ 
tions  are  made  because  the  added  mass  and  moment  of  in¬ 
ertia  of  added-mass  coefficients  are  mainly  dependent  upon 
Che  local  flow  field.  Therefore,  a  two-dimensional  result 


may  be  expected  to  provide  some  qualitative  check  because 
the  ship  is  tong  compared  with  the  length  scale  of  the  mid¬ 
ship  cross-section  in  our  problem  at  hand.  In  nondimen- 
sionalizing  our  two-dimensional  hydrodynamic  coefficients, 
we  assume  that  the  two-dimensional  problem  has  a  unit 
length  so  that  the  nondimensionaluation  of  the 
hydrodynamic  coefficients  will  be  consistent  with  the 
definition  given  in  Equations  (23). 

As  expected,  comparisons  of  pyy  and  P44  in  Figure  3 
and  pj 4  in  Figure  $.  with  pjj,  P44,  and  pM  in  Figure  7, 
show  remarkable  similarities  between  them,  except  that  the 
locations  of  the  spikes  in  two-  and  three-dimensions  are 
somewhat  different:  the  spikes  in  the  values  of  pjj,  P44, 
and  pja  occur  at  vB  «  0.784  in  two-dimensions  and  at  vB 
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Figure  5  —  Coupling  Hydrodynamic  Coefficients,  pjj  »nd  versus  rB,  for  l/B  -  0 


-  0.9  through  1.0  in  three -dimensions.  Comparison  of  the 
have  added -mass  coefficient  in  Figure  3  with  that  in 
Figure  7  also  shows  a  similarity  in  their  qualitative  behavior 
when  vB  >0.5.  However,  when  vB  <0.5,  the  three- 
dimensionality  effect  seems  to  be  more  significant,  because 
the  value  of  pjj  increases  in  three-dimensions,  while  it 
decreases  more  rapidly  in  two-dimensions,  as  vB  decreases 
from  0.5  to  0.2. 

A  close  examination  of  the  local  flow  field  at  the  mid- 
ihip  cross-section  corresponding  to  both  positive  and 
negative  peaks,  for  example,  in  the  sway  added  mass  coef¬ 
ficient  jijj  in  Figures  3  and  7,  shows  that  the  relative  fluid 
motion  under  the  ship  is  out  of  phase  for  the  positive  and 
in  phase  for  the  negative  values  of  added  mass.  In  order  to 
have  a  close  examination  on  the  location  of  these  spikes,  in 
addition,  we  also  treated  an  eigenvalue  problem  correspon¬ 
ding  to  the  foregoing  two-dimensional  boundary-value  pro¬ 
blem  by  a  finite  element  method.  When  the  finite  element 
method  is  applied  to  the  eigenvalue  problem,  the  (original) 
eigenvalue  problem  defined  in  a  partial  differential  equa¬ 
tion  reduces  to  a  general  matrix  eigenvalue  problem  of  a 
type  [A]  X  -  1(C)  X,  where  (A)  and  (Cl  are  matrices  and 
X  is  an  eigenvector.  In  the  present  computations  two  sets 
of  eigenvalue  problems  are  solved.  In  the  first  set,  (A]  and 
|C)  are  identical  to  the  matrices  used  in  solving  the 
boundary-value  problem  to  obtain  the  results  given  in 
Figure  7.  To  obtain  this  set  of  matrices,  96  quadrilateral 
elements  with  a  total  of  345  nodes  in  the  fluid  and  18 
nodes  on  the  free  surface  are  used.  As  a  result,  18  eigen¬ 
values  and  eigenvectors  are  obtained  in  the  first  set.  In  the 
second  set,  the  submatrices  are  taken  for  (A)  and  (C|  from 
the  matrices  constructed  for  our  three-dimensional 


boundary-value  problem  at  hand  by  retaining  the  matrix 
elements  corresponding  to  the  nodes  at  the  midship  cross- 
sectional  plane.  In  the  second  set,  in  which  a  coarse  mesh 
subdivision  (total  of  25  nodes)  is  used,  four  eigenvalues 
and  eigenvectors  are  obtained. 

The  lowest  nonzero  eigenvalue  is  vB  *  0.78410  in  the 
first  set  which  has  fine  mesh  subdivisions  and  vB  » 

0.82912  in  the  second  coarse  mesh  subdivisions.  The  higher 
modes  of  eigenvalues  we  computed  are  out  of  the  range 
considered  here.  However,  the  location  of  the  spikes  occur¬ 
ring  in  our  three-dimensional  results,  for  example  pjj  in 
Figure  3,  is  between  vB  *=  0.9  and  1.0  not  at  vB  = 

0.82912.  If  we  had  solved  a  three-dimensional  eigenvalue 
problem  corresponding  to  the  three-dimensional  boundary- 
value  problem  treated  in  the  present  paper,  we  could  have 
obtained  eigenvalues  all  of  which  have  a  nonzero  imaginary 
part.  However,  the  three-dimensional  eigenvalue  problem  is 
not  treated  here.  The  eigenvalues  of  the  two-dimensional 
problem  treated  here  are  all  real.  From  the  comparisons 
among  the  two  lowest  eigenvalues,  vB  «  0.78410  and 
0.82912,  ar.i  the  location  of  the  spikes  in  the  present  three- 
dimensional  numerical  results  (pjj,  M44,  and  pj4,  i.e.,  vB  is 
between  0.9  and  1 .0),  we  can  conclude  that  the  difference 
in  the  locations  of  the  spikes  in  pj}  in  Figures  3  and  7  is 
caused  partially  by  the  inaccuracy  in  our  three-dimensional 
computations  (due  to  the  use  of  not -fine-enough  mesh  sub¬ 
divisions)  and  partially  by  the  nature  of  three- 
dimensionality. 

It  should  be  kept  in  mind  that  in  the  values  of  p2-»  in 
Figure  7,  If  we  had  computed  p«  •*  »B  -  0.7841.  then  we 
could  have  noticed  a  singular  behavior  like  the  spikes 


We  anticipated  some  singular  behavior  near  vB  « 

1.44,  which  is  the  first  mode  of  the  canal  resonance.  So  we 
made  additional  computations  around  this  value  (excluding 
this  value  as  discussed  earlier)  by  taking  finer  intervals  in 
»B.  Our  computed  results  do  not  show  any  abnormality 
like  spikes,  contrary  to  our  anticipation. 

In  Figures  S  through  10,  the  hydrodynamic  coeffi¬ 
cients,  py  and  ljj.  are  shown  for  the  case  of  an  off-center 
location  of  a  ship  in  a  canal.  s/B  =  0.125.  The  excitation 
forces  and  moments  for  this  case  are  shown  in  Figure  II. 
The  general  behavior  of  all  results  for  an  off-center  case  is 
similar  to  those  previously  shown,  except  that  the  ap¬ 
pearance  of  more  sudden  spikes  is  observed.  The  values  of 
**1 1»  £35«  •od  Mm  in  Figure  8  show  spikes  at  vB  =  1.2  and 
1.3.  The  values  of  pji»  Pjj.  and  H44  have  spikes  at  vB  — 


0.0  0.5  1.0  1.5  2.0  2.5 

vB 

Figure  I  —  Added  Masses  and  Moments  of  Inertia  of 


0.9  and  1.0  in  the  same  figure.  Of  interest  is  that  the  valuta 
Of  Mil.  F22-  “d  H66 in  Figure  8  and  p |g.  1|6.  and  in 
Figure  10  show  spikes  which  were  not  shown  in  the 
previous  case  (s/B  =  0).  This  can  be  interpreted  in  the 
sense  that  the  range  of  the  influence  of  the  local  eigen¬ 
values  is  larger  in  the  case  of  s/B  -  0.125  than  in  the  case 
of  s/B  *  0. 

Similar  spikes  are  also  observed  in  the  damping  coeffi¬ 
cients  in  Figure  9.  The  damping  coefficients  and  I44 
reach  a  maximum  peak  at  vB  =  1.0  and  for  i,j  at  vB  » 
1.2.  Some  of  the  coupling  hydrodynamic  coefficients  are 
shown  in  Figure  10.  The  magnitudes  of  wave  excitation 
forces  and  moments  are  shown  in  Figure  II.  The  behavior 
of  wave  excitation  forces  and  moments  are  very  similar  to 
the  damping  coefficients  corresponding  to  motion  which 
has  already  been  discussed  in  the  previous  case  of  s/B  »  0. 
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0.05 


Figure  10  —  Some  of  Coupling  Hydrodynamic  Coefficients  versus  vB  for  s/B  -  0.12$ 


method,  the  main  core  memory  space  can  be  reduced  by 
half  roughly.  The  Central  Processor  Unit  Time  in  the  CDC 
MOO  was  approximately  40  seconds  to  compute  six  motions 
and  one  diffraction  problem  for  one  wave  number  by  solv¬ 
ing  the  complex  matrix  directly.  By  the  second  method,  ap¬ 
proximately  20  percent  of  the  total  CPU  time  was  saved.  In 
the  present  computer  program,  a  considerable  amount  of 
time  is  used  for  the  input-output  operation  due  to  out  of - 
core  storage  use. 
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